module seed

    use global
    use track, only: track_particle

    implicit none
    private

    public :: seed_horiz, seed_zonal, seed_meridional

    contains

    subroutine seed_horiz(i0, j0, k0, l0, rz0)

        integer, intent(in) :: i0, j0, k0, l0
        real, intent(in) :: rz0
        
        integer  :: i, j, k, l
        real :: rx, ry, rz
        real :: vtot, nn
        integer :: n, ii, jj

        i = i0
        j = j0
        k = k0
        l = l0
        rz = rz0

        ! Total transport at this point.
        vtot = -bz(i,j,k,l)*rz -dz(i,j,k,l)
        ! Number of particles to release, each corresponding to ``transport``
        ! m3/s.
        nn = abs(vtot/transport)
        n = int(sqrt(nn))

        do ii=1,n
        do jj=1,n
            rx = real(ii)/real(n+1)
            ry = real(jj)/real(n+1)
            call track_particle(i, j, k, l, rx, ry, rz)
        enddo
        enddo

    end subroutine seed_horiz

    subroutine seed_zonal(i0, j0, k0, l0, ry0)

        integer, intent(in) :: i0, j0, k0, l0
        real, intent(in) :: ry0

        integer :: i, j, k, l
        real :: rx, ry, rz
        real :: vtot, nn
        integer :: n, ii, kk

        i = i0
        j = j0
        k = k0
        l = l0
        ry = ry0

        ! Total transport at this point.
        vtot = -by(i,j,k,l)*ry -dy(i,j,k,l)
        ! Number of particles to release, each corresponding to ``transport``
        ! m3/s.
        nn = abs(vtot/transport)
        n = int(sqrt(nn))

        do ii=1,n
        do kk=1,n
            rx = real(ii)/real(n+1)
            rz = real(kk)/real(n+1)
            call track_particle(i, j, k, l, rx, ry, rz)
        enddo
        enddo

    end subroutine seed_zonal

    subroutine seed_meridional(i0, j0, k0, l0, rx0)

        integer, intent(in) :: i0, j0, k0, l0
        real, intent(in) :: rx0

        integer :: i, j, k, l
        real :: rx, ry, rz
        real :: vtot, nn
        integer :: n, jj, kk

        i = i0
        j = j0
        k = k0
        l = l0
        rx = rx0

        ! Total transport at this point.
        vtot = -bx(i,j,k,l)*rx -dx(i,j,k,l)
        ! Number of particles to release, each corresponding to ``transport``
        ! m3/s.
        nn = abs(vtot/transport)
        n = int(sqrt(nn))

        do jj=1,n
        do kk=1,n
            ry = real(jj)/real(n+1)
            rz = real(kk)/real(n+1)
            call track_particle(i, j, k, l, rx, ry, rz)
        enddo
        enddo

    end subroutine seed_meridional

end module seed
